{"cells":[{"cell_type":"markdown","metadata":{"collapsed":false,"id":"B2E13C03C0E748919DA40A0B58CB4468","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["# Lecture 10: linear model and model diagnostics \n","\n","## Instructor： 胡传鹏（博士）[Dr. Hu Chuan-Peng]\n","\n","### 南京师范大学心理学院[School of Psychology, Nanjing Normal University]"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"543A644E8F7944B9BAFCD281B8818FDC","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["##  Recap: Linear Model and model diagnositcs\n","\n","* Workflow\n","* MCMC diagnostics"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"4F862DA11B554DB289B4E89D1E2922C4","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["我们在上节课中采用一个简化的workflow，具体包括如下几个步骤：\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D1C645DC1627451E809A160B16451896","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["\n","![Image Name](https://cdn.kesci.com/upload/image/rkvikqg9q6.png?imageView2/0/w/960/h/960)\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"C0D7787982734EC384AD5C3B652E8A0D","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["现在，我们通过上节课的例子来简单回顾一下workflow"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"A2C3D390E20D44858DF0526DD185459C","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["在本例中，涉及：研究问题、数据收集、选择模型、选择先验、模型拟合、采样过程评估、模型诊断"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"6FF8EC3420AA40608BCFA62EF53D10CB","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (1) 提出研究问题\n","\n","还有研究发现个体创新行为可能与自尊水平有关。\n","SES_t:自尊水平；\n","EIB_t创新行为\n","\n","试探究两者之间的关系。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"67488C0B681C40AF8373C9C410AEF7F8","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (2) 数据收集"]},{"cell_type":"code","execution_count":1,"metadata":{"collapsed":false,"id":"8DFB687F35914AFB8C05853D55E5D139","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["WARNING (theano.link.c.cmodule): install mkl with `conda install mkl-service`: No module named 'mkl'\n"]},{"data":{"text/plain":["Text(0, 0.5, 'EIB_t')"]},"execution_count":1,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/8DFB687F35914AFB8C05853D55E5D139/rkyxj65t83.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["#加载需要使用的库\n","%matplotlib inline\n","import numpy as np \n","from scipy import stats\n","import matplotlib.pyplot as plt\n","import pandas as pd\n","import arviz as az\n","import pymc3 as pm\n","# mpl_toolkits.mplot3d是用于三维画图的，Axes3D用于画三维图\n","from mpl_toolkits.mplot3d import Axes3D\n","#数据预处理与可视化\n","np.random.seed(123)  # 随机数种子，确保随后生成的随机数相同\n","data = pd.read_csv(\"/home/mw/input/data9464/clean.csv\") # 读取数据，该数据需要提前挂载\n","data['SES_t'] = (data['SES_t'] - data['SES_t'].mean()) / data['SES_t'].std()#将变量进行标准化\n","\n","data['EIB_t'] = (data['EIB_t'] - data['EIB_t'].mean()) / data['EIB_t'].std()#将变量进行标准化\n","\n","plt.scatter(data['SES_t'],data['EIB_t'])\n","plt.xlabel('SES_t')\n","plt.ylabel('EIB_t')"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"61C5CABC8CA34E878D153FBA69542021","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (3) 选择模型"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"59B8B93E4D324F5DAB64DA0CEB5CCA31","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["尝试构建一个简单的线性模型\n","\n","线性模型可以用概率的形式进行表达\n","\n","\n","$\\alpha \\sim Normal(0,1)$      ->      a $\\sim$ Normal(mu,sigma)\n","\n","$\\beta  \\sim Normal(0,1)$      ->      b $\\sim$ Normal(mu,sigma)\n","\n","$\\sigma \\sim HalfNormal(1)$      ->      sigma $\\sim$ HalfNormal(1)\n","\n","$\\mu_i  = \\alpha + \\beta *x$      ->      mu = alpha + beta*x \n","\n","$y \\sim Normal(\\mu_i,sigma)$      ->      y $\\sim$ Normal(mu,sigma)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"C24612EBE8DE46538C3C68A040CD7B57","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (4)选择先验"]},{"cell_type":"code","execution_count":2,"metadata":{"collapsed":false,"id":"8E89C676474E4E2AA6E63B19D7FA8B90","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["# 在pymc3中，pm.Model()定义了一个新的模型对象，这个对象是模型中随机变量的容器\n","# 在python中，容器是一种数据结构，是用来管理特殊数据的对象\n","# with语句定义了一个上下文管理器，以 linear_model 作为上下文，在这个上下文中定义的变量都被添加到这个模型\n","with pm.Model() as linear_model:\n","    # 先验分布: alpha, beta, sigma这三个参数是随机变量\n","    alpha = pm.Normal('alpha',mu=0,sd=1)\n","    beta = pm.Normal('beta',mu=0,sd=1, shape=1)  \n","    sigma = pm.HalfNormal('sigma',sd=1)\n","    \n","    # 先验预测检查\n","    prior_checks = pm.sample_prior_predictive(samples=50)"]},{"cell_type":"code","execution_count":3,"metadata":{"collapsed":false,"id":"47053D46ABF54C6FB352FD0DF7715D11","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/47053D46ABF54C6FB352FD0DF7715D11/rkyxjpx9zt.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["fig, ax = plt.subplots()\n","x1 = np.linspace(-3, 3, 50) #生成从-2，2之间的50个假数据\n","\n","for a, b in zip(prior_checks[\"alpha\"], prior_checks[\"beta\"]):\n","    y = a + b * x1 # 基于假数据生成预测值\n","    ax.plot(x1,y)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"9DA3B2B9238B492D85420F56F0A3C620","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (5) 拟合数据"]},{"cell_type":"code","execution_count":4,"metadata":{"collapsed":false,"id":"8E7C0A6A910546508B4063B77A8C8B8B","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["linear_model = pm.Model()\n","with linear_model :\n","    # 在pymc3中，pm.Model()定义了一个新的模型对象，这个对象是模型中随机变量的容器\n","    # 在python中，容器是一种数据结构，是用来管理特殊数据的对象\n","    # with语句定义了一个上下文管理器，以 linear_model 作为上下文，在这个上下文中定义的变量都被添加到这个模型\n","    alpha = pm.Normal('alpha',mu=0,sd=1)\n","    beta = pm.Normal('beta',mu=0,sd=1,shape=1)\n","    sigma = pm.HalfNormal('sigma',sd=1)\n","    # x为自变量，是之前已经载入的数据\n","    x = pm.Data(\"x\", data['SES_t'])\n","    # 线性模型：mu是确定性随机变量，这个变量的值完全由右端值确定\n","    mu = pm.Deterministic(\"mu\", alpha + beta*x) \n","    # Y的观测值，这是一个特殊的观测随机变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n","    y_obs = pm.Normal('y_obs',mu=mu,sd=sigma,observed=data['EIB_t'])"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"9DE229BD2025423B85C917275D6C3E63","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (6)采样过程诊断\n","\n","如果使用MCMC对后验进行近似，则需要首先对MCMC过程进行评估。\n","\n","* 是否收敛；\n","* 是否接近真实的后验。\n","\n","对采样过程的评估我们会采用目视检查或rhat这个指标"]},{"cell_type":"code","execution_count":5,"metadata":{"collapsed":false,"id":"7776FCCFB2C5441893BE0B223D3C3B89","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Multiprocess sampling (2 chains in 2 jobs)\n","NUTS: [sigma, beta, alpha]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"name":"stderr","output_type":"stream","text":["Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 6 seconds.\n"]}],"source":["with linear_model :\n","    # 使用mcmc方法进行采样，draws为采样次数，tune为调整采样策略的次数，这些次数将在采样结束后被丢弃，\n","    # target_accept为接受率， return_inferencedata=True为该函数返回的对象是arviz.InnferenceData对象\n","    # chains为我们采样的链数，cores为我们的调用的cpu数，多个链可以在多个cpu中并行计算，我们在和鲸中调用的cpu数为2\n","    trace = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)"]},{"cell_type":"code","execution_count":6,"metadata":{"collapsed":false,"id":"215CDFEF736E45A68BCFFCB03E74A7B6","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["array([[<AxesSubplot:title={'center':'alpha'}>,\n","        <AxesSubplot:title={'center':'alpha'}>],\n","       [<AxesSubplot:title={'center':'beta'}>,\n","        <AxesSubplot:title={'center':'beta'}>],\n","       [<AxesSubplot:title={'center':'sigma'}>,\n","        <AxesSubplot:title={'center':'sigma'}>]], dtype=object)"]},"execution_count":6,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/215CDFEF736E45A68BCFFCB03E74A7B6/rkyxldpabt.png\">"],"text/plain":["<Figure size 864x432 with 6 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["az.plot_trace(trace,var_names=['alpha','beta','sigma'])"]},{"cell_type":"code","execution_count":7,"metadata":{"collapsed":false,"id":"9C6692072E1C44D78988AF65F5E5A4ED","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>mcse_mean</th>\n","      <th>mcse_sd</th>\n","      <th>ess_bulk</th>\n","      <th>ess_tail</th>\n","      <th>r_hat</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>alpha</th>\n","      <td>0.0</td>\n","      <td>0.0</td>\n","      <td>5018.0</td>\n","      <td>2745.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>beta[0]</th>\n","      <td>0.0</td>\n","      <td>0.0</td>\n","      <td>4503.0</td>\n","      <td>3309.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>sigma</th>\n","      <td>0.0</td>\n","      <td>0.0</td>\n","      <td>4367.0</td>\n","      <td>3085.0</td>\n","      <td>1.0</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["         mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat\n","alpha          0.0      0.0    5018.0    2745.0    1.0\n","beta[0]        0.0      0.0    4503.0    3309.0    1.0\n","sigma          0.0      0.0    4367.0    3085.0    1.0"]},"execution_count":7,"metadata":{},"output_type":"execute_result"}],"source":["az.summary(trace, var_names=['alpha','beta','sigma'], kind=\"diagnostics\")"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"0C5F6A4325D84CF58C46E8E009EB0CA3","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### (7)模型诊断\n","\n","在MCMC有效的前提下，需要继续检验模型是否能够较好地拟合数据。\n","\n","我们会使用后验预测分布通过我们得到的参数生成一批模拟数据，并将其与真实数据进行对比。"]},{"cell_type":"code","execution_count":8,"metadata":{"collapsed":false,"id":"08DBC6B10B134772855BE0F0CF010E27","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"data":{"text/html":["\n","    <div>\n","      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n","      100.00% [4000/4000 00:05&lt;00:00]\n","    </div>\n","    "],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"update_display_data"},{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n","  FutureWarning,\n"]}],"source":["# 后验预测分布的计算仍在容器中进行\n","with linear_model:\n","    # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n","    ppc_y = pm.sample_posterior_predictive(trace.posterior) \n","#将ppc_y转化为InferenceData对象合并到trace中\n","az.concat(trace, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)"]},{"cell_type":"code","execution_count":9,"metadata":{"collapsed":false,"id":"B3201FD32B234C05BD6942227CDE422E","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["<AxesSubplot:xlabel='y_obs'>"]},"execution_count":9,"metadata":{},"output_type":"execute_result"},{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/IPython/core/events.py:89: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n","  func(*args, **kwargs)\n"]},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/B3201FD32B234C05BD6942227CDE422E/rkyxmeilvx.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# 绘制后验预测分布\n","az.plot_ppc(trace)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"57BE24124330453D9B958A8CA14CDB48","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["## Part 2: Model comparation"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"9BF407CFF89149FABD2C42487A5E9ADA","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["我们在上节课中提到的workflow具体包括如下几个步骤：\n","\n","研究问题、选择模型、选择先验、模型拟合、采样过程评估、模型诊断、模型比较、统计推断、结果报告。\n","\n","![Image Name](https://cdn.kesci.com/upload/image/rkm3pw954u.png?imageView2/0/w/960/h/960)\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D69BCF4358754FB286D2FEA8C4A6F149","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["在贝叶斯的workflow中，模型诊断就是对MCMC进行检验，检验的核心在于确定MCMC算法采样得到的样本是否足以提供目标分布的精确近似。\n","\n","模型诊断结果表明MCMC算法采样得到的样本能够提供目标分布的精确近似之后，我们需要确定模型既不太简单以至于错过数据中有价值的信息(underfitting)，也不会太复杂从而过拟合数据中的噪音(overfitting)。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"B549FBCA33B34ECAB8A4403241264779","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["\n","![Image Name](https://img-blog.csdnimg.cn/20210303141100499.png#pic_center)\n","\n","资料来源：https://blog.csdn.net/weixin_43378396/article/details/90707493"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"7C33039E93CF4261A664628B5D5B9D68","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### 欠拟合(underfitting)\n","\n","欠拟合的模型在当前样本的数据拟合效果不好，且其泛化能力(模型在当前样本外新的数据上的预测的准确度)也同样不佳。\n","\n","### 过拟合(overfitting)\n","\n","过拟合的模型在当前样本的数据上的拟合程度很好，但是泛化能力也较差。\n","\n","这是因为模型把训练样本学习地“太好了”，把样本自身地一些特点也当作了所有潜在样本都会具有的一些性质，这样就会导致其泛化性能下降。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"FD588A1708B8432D857C2BF7F362D2FC","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["\n","### 为什么会发生欠拟合\n","\n","- 数据特征较少\n","\n","数据特征指的是数据的属性，比如part1中展示的数据的各个变量就是数据的特征。在所有变量都能独立地对目标变量做出解释的前提下，数据特征越多，数据拟合程度越好。\n","\n","- 模型复杂度过低\n","\n","模型的复杂度代表模型能够描述的所有函数，比如线性回归最多能表示所有的线性函数。\n","\n","模型的复杂度和模型的参数数量有关，一般来说，模型参数越多，复杂度越高，模型参数越少，复杂度越低。\n","\n","### 为什么会发生过拟合\n","\n","- 当前样本的噪音过大，模型将噪音当作数据本身的特征\n","\n","当数据的有些特征与目标变量无关，这些特征就是噪音，但它也可能被误当作数据特征，这就会造成模型过拟合\n","\n","- 样本选取有误，样本不能代表整体\n","\n","- 模型参数太多，模型复杂度太高\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"442CF350FF804ACE8048EEAA95234CF8","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["\n","![Image Name](https://vitalflux.com/wp-content/uploads/2020/12/overfitting-and-underfitting-wrt-model-error-vs-complexity-768x443.png)\n","\n","资料来源：https://vitalflux.com/overfitting-underfitting-concepts-interview-questions/"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"6183013BDE0546EBA62A7B6501434279","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### 如何避免欠拟合\n","\n","- 增加数据的特征\n","\n","- 增加模型复杂度\n","\n","### 如何避免过拟合\n","\n","- 选择更具代表性的数据\n","\n","- 降低模型复杂度"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"942ACCF889A549AA8A3E55D0283DEC0B","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["为了选择一个能够在过拟合和欠拟合之间的达到平衡的最佳模型，就需要进行模型比较。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"26E115511E08465DABDC1532645E74F0","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["根据不同模型比较的指标的关注点和原理，进行模型比较的统计指标可被分为三类。 \n","\n","第一类是模型的拟合优度(Goodness of fit)\n","\n","第二类是交叉验证(Cross validation)和近似交叉验证\n","\n","第三类是边缘似然理论的指标"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"A0EB6C8163C44F2BA5298F031068F84D","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["\n","\n","### 拟合优度(Goodness of fit)\n","\n","拟合优度衡量模型在样本集上拟合的程度，并没有考虑模型的复杂度，只是单纯地用于衡量模型在当前样本的拟合的好坏。\n","\n","比如回归模型中的$R^2$，即回归模型所能解释的因变量的变异的百分比。\n","\n","常用的指标有均方差(Mean squared deviation, MSD)，对数似然((log likelihood function)\n","\n","\n","\n","![Image Name](https://th.bing.com/th/id/R.f8319827592edfae31cd1248ebb7be01?rik=UksWyJLwHYJqXQ&pid=ImgRaw&r=0)\n","\n","\n","\n","资料来源：https://baike.baidu.com/pic/%E6%8B%9F%E5%90%88%E4%BC%98%E5%BA%A6/10946257/0/6a63f6246b600c338b38d998114c510fd9f9a118"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"F142E3885204493B853869C751D9FD4D","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 均方差\n","\n","计算公式为：$MSE(\\hat{\\theta})=E(\\hat{y}-y)^2=\\frac{1}{n}\\sum (\\hat{y}-y)^2$\n","\n","均方差表示预测变量与目标变量之间差值的平方之和的均值\n","\n","\n","![Image Name](https://cdn-media-1.freecodecamp.org/images/MNskFmGPKuQfMLdmpkT-X7-8w2cJXulP3683)\n","\n","资料来源：https://www.freecodecamp.org/news/machine-learning-mean-squared-error-regression-line-c7dde9a26b93/\n","\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"BD33902DD7F44B3B963982F2EC3B70B6","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["MSE与许多常见的拟合优度指标有关联，\n","\n","对MSE开根号，可以得到**均方根误差(Root mean square deviation, RMSD)**;\n","\n","给MSE乘上数据点的数量，即可得到**残差平方和(Residual sum of squares, RSS)** 。\n","\n","MSE在认知建模中往往应用在对连续变量的数据的建模上，另外RSS还可以在计算AIC和BIC时替代对数似然函数。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"03E73CE475354F568C7282529BA1AFA8","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 对数似然(Log-likelihood)\n","\n","\n","\n","对数损失似然函数的公式为：\n","\n","$L=\\prod (P_i) = p(x_1) * p(x_2) * ...* p(x_n)$ \n","\n","似然(likelihood)表示每个数据在当前模型下发生的可能性的乘积\n","\n","$LL =log(L)=\\sum log(P_i)$ \n","\n","由于乘积的形式会得到一个很小的数，不便于后续处理，因此对数似然将似然取对数以便比较。\n","\n","\n","横轴为待评估的模型的各参数值，左侧纵轴为似然的值，右侧纵轴为对数似然的值\n","\n","\n","![Image Name](https://www.aptech.com/wp-content/uploads/2020/09/poisson-likelihood-function.jpeg)\n","\n","资料来源：https://www.aptech.com/blog/beginners-guide-to-maximum-likelihood-estimation-in-gauss/\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"8D943CA781014C7D91628862D3AE5647","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["对数损失似然函数有一大优点，即在不论建模的数据为离散分布的选项数据还是连续分布的反应时或者评分数据时等情况都可以使用。\n","\n","在计算模型的研究里，对数似然函数作为模型拟合优度的指标有两种用途。\n","\n","第一，研究者可以对比不同模型的平均对数似然函数来判断模型是否对数据有很好的拟合。\n","\n","第二，研究者可以使用似然比检验(Likelihood test)去判断嵌套模型之间的差异是否显著。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"AB2F72934DA24227B3A2E1EF18DAD1C4","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### 模型的泛化能力：预测准确性\n","\n","前面提到的诸多模型拟合优度 (Goodness of fit) 指标，只能衡量模型对于当前样本的拟合程度。\n","\n","对于样本外的数据，我们不确定该模型是否具有**泛化能力**，即该模型是否能准确的预测样本外的数据。\n","\n","为了评估模型的预测能力，我们有以下3种策略：\n","1. 通过新数据对模型进行评估。我们可以收集新的数据，并检验模型的预测能力。\n","2. 从已有样本中拿出一部分数据用来预测。即交叉验证。\n","3. 通过统计方法进行近似。比如使用对交叉验证近似的信息熵指标来评估模型的泛化能力。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"449A531CF85F4E9ABEED728BB22EABF7","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 交叉验证\n","\n","收集新数据来检验模型的预测能力是一种理所当然的直觉。但心理学数据不同于其他学科的数据，它常受到**时间因素**和**抽样**的影响。\n","- 比如，心境可能随着季节变化，因此在不同季节收集到的数据会受到时间的影响。\n","\n","因此，一种更高效的方法是，一次性多收集一些数据，选择其中的一部分作为预测数据。\n","\n","但问题在于，我们选择哪一部分数据作为预测数据呐？或者说，我们该如何有效的对数据进行抽取呐？\n","\n","**交叉验证**的目的就在于：提供不同的抽取预测数据的策略"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"B1CADF6BDCDC42E79918C6A1FC066CBC","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["常见的交叉验证策略：\n","1. 分半交叉验证 (Split-half cross-validation)\n","\t- 分半交叉验证将观测数据对半分成两部分，分别在不同的数据集上拟合模型，并在另外一半数据集上验证模型，最后再对比不同的模型在两份数据集作为验证集时的预测准确度。\n","2. K 折交叉验证 (K-fold cross-validation) \n","\t- K 折交叉验证把数据分成 K 分，其中一份作为训练集，其余的 K-1 分数据集作为验证集，总共重复这个流程 K 次。以 K 次验证结果的均值作为验证标准。\n","3. 留一法交叉验证 (Leave-one-out cross-validation)\n","\t- 留一法交叉验证是 K 折交叉验证的一个特例，当分折的数量等于数据的数量时，K 折留一法便成了留一法交叉验证。留一法交叉验证相较于普通的交叉验证方法，几乎使用了所有数据去训练模型，因此留一法交叉验证的训练模型时的**偏差 (bias) 更小、更鲁棒**，但是又因为验证集只有一个数据点，验证模型的时候**留一法交叉验证的方差 (Variance) 也会更大**。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"268B2411BC6A473E9DE411265BF42AAC","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["K 折交叉验证图示\n","![](https://hub.packtpub.com/wp-content/uploads/2019/05/KFold.png)\n","\n","source: https://hub.packtpub.com/cross-validation-strategies-for-time-series-forecasting-tutorial/"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"15E81112F633419FA4F9DB96C71AC219","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 信息准则(information criteria)\n","\n","尽管交叉验证存在诸多优势，并且避免了收集数据的潜在问题，但是交叉验证**在认知神经科学里的认知建模领域里的使用并不广泛**，\n","\n","其最主要的原因在于，研究者收集的数据**样本量往往有限**，然而很多计算模型却对数据样本量有所需求，如果使用交叉验证，将数据拆分，那么很有可能导致拟合模型的试次数量不足，使得模型拟合和验证的结果较差，进而导致产生一类错误和二类错误的概率增大。\n","\n","为了解决这些问题，统计学家们凭借数学与统计工具可以近似得到与交叉验证相似的结果。\n","\n","这类近似交叉验证的指标被称为**信息准则**(information criteria)。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"A835429E78A748B5BE97549F139A2844","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["信息准则是对统计模型**预测精度**的一种度量。它考虑**模型与数据的匹配程度**，并通过**模型的复杂性**进行惩罚。\n","\n","从前面关于模型**偏差方差权衡** (Bias variance trade-off)的讨论中我们了解到\n","- **偏差 bias** 会随着模型的复杂度的增大而减小，偏差大的模型往往*欠拟合*；\n","- **方差 error** 则会随着模型复杂度的增大而增大，方差大的模型往往*过拟合*。\n","- 方差和偏差之间存在一种**权衡关系**。\n","- 一般来说，复杂的模型容易过拟合，而简单的模型则容易欠拟合。\n","\n","其中：\n","- **模型与数据的匹配程度**更多的体现在模型的**偏差 bias**上。\n","- 而**模型的复杂性**更多的体现在模型的**方差 error**上。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"A3BA6B8A87B3496C9D44EF22D9BA5588","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["通过稍微数学化的语言表达就是：\n","\n","$信息准则 = 偏差 + 方差 = bias + error$\n","- 其中 bias 往往与模型的似然值相关。可以想象，似然值越大，模型的偏差越小。\n","- error 与模型的复杂程度相关。比如，模型的参数越多模型越复杂。\n","\n","可见，信息准则越小，偏差和方差就越小。因此，**信息准则越小，代表模型的预测性越好**。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"F22BA4CC600F4E28B125358CAD77ACB1","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["常见的4类信息准则：\n","1. AIC (Akaike information criterion)\n","2. DIC (Deviance information criterion)\n","3. WAIC (Widely applicable information criterion)\n","4. LOO-CV (Leave-one-out cross-validation)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"63741BADA01A47B294A0839838677898","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 1. AIC (Akaike information criterion)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"7ABEB248AF6F4249831570CE2472A8D4","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["AIC是最简单的信息准则标准，由日本统计学家赤池弘次 (Hirotugu Akaike) 提出 (Akaike, 1974)，是频率主义统计学里最为经典的模型比较指标之一。\n","\n","其表达式如下：\n","$$\n","A I C= bias + error = -2 \\sum_{i}^{n} \\log p\\left(y_{i} \\mid \\hat{\\theta}_{m l e}\\right)+2 p_{A I C}\n","$$\n","\n","别被公式吓到🧐，它的原理其实很简单。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"7FB883B7C30141078092183CC90662E5","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["$A I C= bias + 2 p_{A I C}$\n","\n","- 其中，$2 p_{A I C}$ 就是 error，代表模型参数数量的两倍。\n","- error 描述了模型的复杂程度。即参数越多，模型越复杂，相应的 AIC 值越大，模型预测性越差。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"316727DBDAB14FB6988BEA5AFE02591C","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["$A I C = -2 \\sum_{i}^{n} \\log p\\left(y_{i} \\mid \\hat{\\theta}_{m l e}\\right)+2 p_{A I C}$\n","\n","- $-2 \\sum_{i}^{n} \\log p\\left(y_{i} \\mid \\hat{\\theta}_{m l e}\\right)$ 为 bias，描述了模型对于当前数据的匹配程度。\n","- bias 与 似然函数 $p\\left(y_{i} \\mid \\hat{\\theta}_{m l e}\\right)$ 相关。$\\hat{\\theta}_{m l e}$为最大似然法求得的参数值。 \n","  \n","  我们知道，模型拟合的越好，似然值越大，因此$\\sum_{i}^{n} \\log p\\left(y_{i} \\mid \\hat{\\theta}_{m l e}\\right)$越大。相应的 $-2\\sum_{i}^{n} \\log p\\left(y_{i} \\mid \\hat{\\theta}_{m l e}\\right)$ 的值越小。那么 AIC的值也越小。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"1969289189F7477D81C4D7AAB053F640","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["需要注意的是 AIC 只考虑了最大似然对应的参数值 $\\hat{\\theta}_{m l e}$, 因此它适用于频率学派模型的评估。而对于贝叶斯学派来说，**由于参数为分布，因此不能使用AIC来评估模型**。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"F6A3B2A44A7E436A82E7203F999F3247","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 2. DIC (Deviance information criterion)\n","\n","为了解决 AIC 无法评估贝叶斯模型。加之 AIC 只考虑到了**模型参数数量**所带来的复杂性。\n","\n","统计学家们提出了\"贝叶斯参数估计版的 AIC\"，即 DIC (Deviance information criterion) 。\n","\n","$$\n","\\mathrm{DIC} = -2 \\sum_{i}^{n} \\log p\\left(y_{i} \\mid \\bar{\\theta}\\right) +2 p_{D}\n","$$"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"1CF12B426CF1487B911DDC8733679EE6","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["可以看到，DIC 与 AIC 的公式非常相似。\n","\n","- 其中$-2 \\sum_{i}^{n} \\log p\\left(y_{i} \\mid \\bar{\\theta}\\right)$ 为 bias。也称为 deviance，即 DIC 的 D。\n","- 与 AIC 的区别在于，由于贝叶斯模型中的参数 $\\theta$ 为概率分布，DIC 选择用参数分布的均值 $\\bar{\\theta}$ 来替代 AIC 中的 $\\hat{\\theta}_{m l e}$。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"6AD7CD1D11794FA6936075E4C9555E8B","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["\n","![Image Name](https://cdn.kesci.com/upload/image/rkvkoshdci.png?imageView2/0/w/640/h/640)\n","\n","DIC 通过 deviance 评估模型的预测性，避免了 AIC 只考虑参数数量复杂性对于模型的影响。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E4D2CC031C994610BD2798D372AC15C1","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["DIC 对于 error 的理解与 AIC 不同。DIC 考虑了参数后验分布中所有采样对于数据的预测性。\n","\n","$$\n","\\frac{1}{2}error = p_{D} = \\bar{D} - D_{\\hat{\\theta}}\n","$$\n","\n","- 要理解 DIC 中的 error 与 DIC 如何评估后验分布中不同参数采样的影响，关键在于了解 D (deviance) 的意义。\n","- D (deviance) 的公式为 $-2 \\sum_{i}^{n} \\log p\\left(y_{i} \\mid \\theta\\right)$。\n","- 可见，bias 也可以用 D 来表示，即 $D_{\\hat{\\theta}}$ 。差异在于，公式中的 $\\theta$ 变为 $\\hat{\\theta}$。 并且这个 bias $\\hat{\\theta}$ 也出现在 $p{D}$ 的公式中。\n","- 有了 $\\hat{\\theta}$，只需要再得到 $\\bar{D}$ 就行了。我们可以用参数后验分布中所有采样(如n个采样)求得，n个 $D_i$， i为1到n。那么这 n个 $D_i$ 的均值就可以表示**参数后验分布中所有采样对于数据的平均预测性**，即 $\\bar{D}$。\n","- $p_{D} = 所有参数采样 D 的均值 - 参数分布均值的D = \\bar{D} - D_{\\hat{\\theta}}$\n","- 可以想象，当两种 D 的差异最小时，$p_{D} =0$，此时模型的预测性最好，DIC值也更低。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"98138646A7D5462690974645D4036660","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["DIC  是心理学领域最常用的模型评估指标之一。\n","- 与 AIC 只是简单地使用了参数数量作为复杂度的惩罚项不同的是，DIC **利用了 MCMC 采样得到的计算模型参数的后验分布去计算模型的有效数量**。\n","- 另外 DIC的计算速度与AIC一样很快，这与后面会介绍的其他指标形成对比。\n","- 最后，除了使用参数分布的均值去计算 DIC 中 bias，也可以**使用参数分布的中位数等计算 bias，这提高了 DIC 计算的灵活性**。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"26F33DA106EC494DA0F7A8ABAEE7A02E","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 3. WAIC (Widely applicable information criterion)\n","\n","DIC 很好的解决了 AIC 无法适用于贝叶斯模型的问题，并且考虑了参数后验分布的影响。\n","\n","DIC 的问题在于：\n","- 如果参数后验分布不是正态分布，那么 $D_{\\hat{\\theta}}$ 的估计会存在偏差。\n","- 这进而反应为，$p{D}$ 是负数。$p{D}$ 为负数导致整个 DIC 为负数，也是 DIC 的一大问题。\n","\n","为了解决上面的问题，日本统计学家渡边澄夫提出了“DIC 的升级版”的升级版，WAIC (Widely applicable information criterion)。\n","\n","$$\n","W A I C=\\sum_{i}^{n} \\log \\left(\\frac{1}{s} \\sum_{j}^{S} p\\left(y_{i} \\mid \\boldsymbol{\\theta}^{j}\\right)\\right)-\\sum_{i}^{n}\\left({V}^{s}_{j} \\log p\\left(Y_{i} \\mid \\boldsymbol{\\theta}^{j}\\right)\\right)\n","$$"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"832D98262DB94193AF874E9F4CA43F2B","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["![Image Name](https://cdn.kesci.com/upload/image/rkvkzpct7z.png?imageView2/0/w/482/h/480)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"826C805F8C7F4C06AFA50144FAE507C3","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["$\\sum_{i}^{n} \\log \\left(\\frac{1}{s} \\sum_{j}^{S} p\\left(y_{i} \\mid \\boldsymbol{\\theta}^{j}\\right)\\right)$ 为 bias 与 似然函数有关。\n","- $\\sum_{i}^{n} \\log \\left(\\frac{1}{s} \\sum_{j}^{S} p\\left(y_{i} \\mid \\boldsymbol{\\theta}^{j}\\right)\\right) = \\sum_{i}^{n} \\log \\bar{L}$，L 为似然函数。\n","- $\\bar{L}$ 类似于 $\\bar{D}$，是参数后验分布中所有参数对应的似然值的均值。\n","- 区别在于：\n","  - WAIC 每次选择一个数据点，计算它在所有后验采样上的似然值 $\\sum_{j}^{S} p\\left(y_{i} \\mid \\boldsymbol{\\theta}^{j}\\right) = \\sum_{i}^{n} \\log \\bar{L}$，\n","  - 再求这些似然值在不同后验上的平均值 $\\frac{1}{s} \\sum_{j}^{S} p\\left(y_{i} \\mid \\boldsymbol{\\theta}^{j}\\right) = \\sum_{i}^{n} \\log \\bar{L}$，\n","  - 最后将不同数据点上的似然值求和，即 $\\sum_{i}^{n} \\log \\left(\\frac{1}{s} \\sum_{j}^{S} p\\left(y_{i} \\mid \\boldsymbol{\\theta}^{j}\\right)\\right) = \\sum_{i}^{n} \\log \\bar{L}$。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"20285AAC5F27430CBC659F9E7878CFD9","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["$\\sum_{i}^{n}\\left({V}^{s}_{j} \\log p\\left(Y_{i} \\mid \\boldsymbol{\\theta}^{j}\\right)\\right)$  为 error，与 bias 不同的是，bias 计算的是似然的均值，而 error 通过 V 计算似然值的变异。\n","- 与上面的流程类似，首先选择一个数据点i，求得它在不同参数j下的似然值 $\\log p(Y_{i} \\mid \\boldsymbol{\\theta}^{j})$，\n","- 计算该数据点在不同参数采样j下似然的变异 $({V}^{s}_{j} \\log p(Y_{i} \\mid \\boldsymbol{\\theta}^{j}))$\n","- 最后，将每个数据点的变异加起来 $\\sum_{i}^{n}({V}^{s}_{j} \\log p(Y_{i} \\mid \\boldsymbol{\\theta}^{j}))$"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"4BB20F8A4C8343D7ACEEFE33934CCB87","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 4. LOO-CV (Leave-one-out cross-validation)\n","\n","WAIC 是非常优秀且常用的指标，其本质是模型对于未知数据的预测能力的**近似**。\n","\n","另一种与 WAIC 非常类似的近似方法是前文提到的**留一交叉验证法** (Leave-one-out cross-validation, LOO-CV) \n","\n","$$\n","ELPD_{LOO-CV} = \\sum_{\\mathrm{y}} \\log E_{\\theta}\\left(L\\left(y_{i}|\\theta_{-i}\\right)\\right)\n","$$"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D040BC3FBBD146BA9E557E5145D48374","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["$ELPD_{LOO-CV}$ 的数学表达式非常类似，其特点为：\n","- $ELPD_{LOO-CV}$ (expected log pointwise predictive density) 利用了留一交叉验证的思想，用去除数据点i剩下的数据拟合模型；再回过来用该模型的参数去预测数据点i。\n","- 因此，ELPD 与 WAIC 的 bias 非常类似。区别在于，似然中的参数值不是通过所有数据进行拟合的，而是通过去除数据点i剩下的数据拟合得到的参数 $\\theta_{-i}$。\n","- 此外，更巧妙的是，由于 $ELPD_{LOO-CV}$ 是直接对未来数据进行预测，因此不需要再惩罚模型，即不需要 error项了。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"4D961BE79CBD470D9D4F433F03968BFB","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["\n","![Image Name](https://cdn.kesci.com/upload/image/rkvl0j94s7.png?imageView2/0/w/640/h/640)\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"E05B1ACC376F44319A541C7CA6225850","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["在实际操作中，我们通过 `ArViz` 的函数可以很容易的获得 WAIC 和 $ELPD_{LOO-CV}$，我们称为 **LOO** 方法。\n","\n","此外，由于 $ELPD_{LOO-CV}$ 的计算量也比较大，ArViz 会使用 Pareto Smooth Importance Sampling Leave Once Out Cross Validation (PSIS-LOO-CV) 来近似 $ELPD_{LOO-CV}$。\n","\n","PSIS-LOO-CV 有两大优势：\n","1. 计算速度快，且结果稳健\n","2. 提供了丰富的模型诊断指标"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"BAC3D301FD7A4A4AB3F901150EE94C98","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["### Model Averaging\n","\n","前面我们讨论了如何评估单个模型的**拟合优度**与**预测精度**。\n","\n","但判断模型预测性能的另一个思路是：拟合多个模型，比较不同模型的预测能力。在实践中，我们对不同的模型赋予不同的权重，并组合他们生成一个元模型 (meta-model)，进而进行元预测，以此评估不同模型权重的影响。\n","\n","这种给不同模型赋予权重的方法称为 模型平均法 Model Averaging。\n","常见有三种计算模型权重的方式：\n","1. Marginal likelihood\n","2. BIC (Bayesian information criterion)\n","3. Pseudo Bayesian model averaging"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"2CE5C18A893341EF88BEE954BE7ED75C","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 1. Marginal likelihood\n","\n","边缘似然或者边际似然 (marginal likelihood) 是贝叶斯公式的分布部分\n","- 即 $p(\\theta|data)=\\frac{p(data|\\theta)p(\\theta)}{p(data)}$ 中的 $p(data) = \\int_{\\theta}^{} p(data|\\theta)p(\\theta)d\\theta$ \n","- 边缘似然与贝叶斯公式分子部分的似然不同，表达了模型对数据的平均拟合 (Average fit)，因此它可以作为模型选择的指标。*边缘似然越大，说明模型对样本数据解释的越好*。\n","- 当比较两个模型时，可以将边缘似然转化为**贝叶斯因子**(Bayes factor)。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"9252096431E94EB8B0709D86806D7304","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["Bayesian Model Averaging\n","\n","$$\n","w_k = \\frac{e^{-ML_{k}}}{\\sum e^{-ML_{k}}}\n","$$\n","\n","通过边缘似然可以计算 Bayesian Model Averaging。\n","- 假设有 k 各模型。\n","- k个模型的边缘似然为 $ML_{k}$。\n","- k个模型的权重 $w_k$ 为当前模型的边缘似然 $ML_{k}$ 比上 所有模型边缘似然之和 $\\sum e^{-ML_{k}}$。\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"FB0E8CFCEA9D41648EA1FBD5E38F844F","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 2. BIC (Bayesian information criterion)\n","\n","因为边缘似然计算量巨大，因此我们需要一些快速的计算方式，比如 BIC (Bayesian information criterion)。\n","\n","$$\n","A I C=-2 \\sum_{i}^{n} \\log p\\left(y_{i} \\mid \\hat{\\theta}_{m l e}\\right)+2 k\n","$$\n","\n","$$\n","B I C=-2 \\sum_{i}^{n} \\log p\\left(y_{i} \\mid \\hat{\\theta}_{m l e}\\right) + 2 k*ln(n)\n","$$\n","\n","BIC 与 AIC 非常类似，区别在于惩罚项 error的不同。\n","- BIC 中的 error 在 AIC error 的基础上乘以 ln(n)。\n","- 其中，k为模型参数的数量；n为数据的数量。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"DFFEF5A074ED4E95B4D5A765CFD44CC9","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["BIC 的特点：\n","- 在公式上与 AIC 高度相似，因此可以用来检验模型拟合优度，其值越小，模型拟合越好。\n","- BIC 的 error 往往比 AIC 的更大，即惩罚更大，因此，**BIC 通常会选择简单的模型**。\n","- BIC 虽然适用于贝叶斯模型，但是它没有考虑先验的影响。\n","- 最重要的是，BIC **是边缘似然的近似**，计算速度比 ML 更快，并且同样也可以被用来计算贝叶斯因子。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"3B6D6D3C66DA4A71834A49A9C31DAEF4","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### 3. Pseudo Bayesian model averaging\n","\n","BIC 虽然能近似边缘似然， 但是其 error 只考虑了参数数量和数据数量的复杂性，也没有考虑到先验的影响。\n","这样计算的模型权重很可能存在偏差。\n","\n","为了更高效的计算模型的权重。一种可行的方法是 Pseudo Bayesian model averaging。\n","- 即通过 WAIC 与 LOO 来近似边缘似然。\n","- 再通过 Bayesian model averaging 公式计算模型权重。\n","\n","\n","$$\n","w_i = \\frac{e^{-\\Delta_{i}}}{\\sum_{j}^{k} e^{-\\Delta_{j}}}\n","$$\n"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"9B858DB4D2AC49628342C875DC189A4D","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["上面的公式与 Bayesian model averaging 的公式 $w_k = \\frac{e^{-ML_{k}}}{\\sum e^{-ML_{k}}}$ 很像。\n","区别在于：\n","- 通过 WAIC 或者 LOO 在模型中的差值 $\\Delta_{j}$ 替代了 边缘似然 ML。\n","- $\\Delta$ 表示的是，第 i 或者 j个模型的 WAIC 与 最优模型的 WAIC 的差值。\n","\n","PyMC3 与 Arviz 提高了很多关于 Pseudo Bayesian model averaging 计算的方法，之后的实践中，我们将注重于通过 Pseudo Bayesian model averaging 展示模型平均法的作用。"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"A65A8202204D49F59C88220EC39413DB","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["总结：\n","\n","模型拟合优度的方法包括：\n","- 拟合优度 \n","- mse \n","- 对数似然\n","\n","模型预测进度的方法包括：\n","- AIC\n","- DIC\n","- WAIC\n","- LOO\n","\n","模型平均法包括：\n","- Bayesian model averaging\n","- BIC\n","- Pseudo Bayesian model averaging"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"DBE02B0FAB6140E383E8301C9FF4EEC6","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["|          | AIC                          | DIC                              | WAIC   | LOOCV  | BIC                          |\n","| -------- | ---------------------------- | -------------------------------- | ------ | ------ | ---------------------------- |\n","| 适用框架 | 频率论                       | 贝叶斯                           | 贝叶斯 | 贝叶斯 | 贝叶斯/频率论                |\n","| bias     | 最大似然所对应参数下的似然值 | 参数均值对应的似然 D               | ELPD   | $ELPD_{LOO-CV}$   | 最大似然所对应参数下的似然值 |\n","| error    | 参数数量的两倍               | $D = \\bar{D} - D_{\\hat{\\theta}}$                      | 似然的变异   |    | 参数数量+数据数量            |"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"73236A5EDD544D3A8C0F1BED5967C88D","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["## Example"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"7688854B622740BE95D71439160FF467","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["认知控制是人的一项重要的能力，在高一致性冲突和低一致性冲突的条件下，反应时间会有不同的表现。\n","\n","我们建立三种不同的模型，来解释探索不同冲突条件下反应时间的分布。\n","\n","数据来源：James F Cavanagh, Thomas V Wiecki, Michael X Cohen, Christina M Figueroa, Johan Samanta, Scott J Sherman, and Michael J Frank. Subthalamic nucleus stimulation reverses mediofrontal influence over decision threshold. Nature Neuroscience, September 2011. URL: http://www.nature.com/doifinder/10.1038/nn.2925, doi:10.1038/nn.2925."]},{"cell_type":"code","execution_count":10,"metadata":{"collapsed":false,"id":"180DED04984E420D8DE2E055D113EBEA","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["import pandas as pd\n","import numpy as np\n","import arviz as az\n","import pymc3 as pm\n"]},{"cell_type":"code","execution_count":14,"metadata":{"collapsed":false,"id":"EB6236DAF4364ED6988BFDA6C6E4731D","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>subj_idx</th>\n","      <th>stim</th>\n","      <th>rt</th>\n","      <th>response</th>\n","      <th>theta</th>\n","      <th>dbs</th>\n","      <th>conf</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>0</th>\n","      <td>0</td>\n","      <td>LL</td>\n","      <td>1.21</td>\n","      <td>1.0</td>\n","      <td>0.656275</td>\n","      <td>1</td>\n","      <td>HC</td>\n","    </tr>\n","    <tr>\n","      <th>1</th>\n","      <td>0</td>\n","      <td>WL</td>\n","      <td>1.63</td>\n","      <td>1.0</td>\n","      <td>-0.327889</td>\n","      <td>1</td>\n","      <td>LC</td>\n","    </tr>\n","    <tr>\n","      <th>2</th>\n","      <td>0</td>\n","      <td>WW</td>\n","      <td>1.03</td>\n","      <td>1.0</td>\n","      <td>-0.480285</td>\n","      <td>1</td>\n","      <td>HC</td>\n","    </tr>\n","    <tr>\n","      <th>3</th>\n","      <td>0</td>\n","      <td>WL</td>\n","      <td>2.77</td>\n","      <td>1.0</td>\n","      <td>1.927427</td>\n","      <td>1</td>\n","      <td>LC</td>\n","    </tr>\n","    <tr>\n","      <th>4</th>\n","      <td>0</td>\n","      <td>WW</td>\n","      <td>1.14</td>\n","      <td>0.0</td>\n","      <td>-0.213236</td>\n","      <td>1</td>\n","      <td>HC</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["   subj_idx stim    rt  response     theta  dbs conf\n","0         0   LL  1.21       1.0  0.656275    1   HC\n","1         0   WL  1.63       1.0 -0.327889    1   LC\n","2         0   WW  1.03       1.0 -0.480285    1   HC\n","3         0   WL  2.77       1.0  1.927427    1   LC\n","4         0   WW  1.14       0.0 -0.213236    1   HC"]},"execution_count":14,"metadata":{},"output_type":"execute_result"}],"source":["data = pd.read_csv('/home/mw/input/data7447/cavanagh_theta_nn.csv') # 载入数据\n","data.head(5) # 查看前5行数据"]},{"cell_type":"code","execution_count":15,"metadata":{"collapsed":false,"id":"F108F2B40F2044458884849B35672501","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>count</th>\n","      <th>mean</th>\n","      <th>std</th>\n","      <th>min</th>\n","      <th>25%</th>\n","      <th>50%</th>\n","      <th>75%</th>\n","      <th>max</th>\n","    </tr>\n","    <tr>\n","      <th>conf</th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","      <th></th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>HC</th>\n","      <td>1972.0</td>\n","      <td>1.454450</td>\n","      <td>0.706737</td>\n","      <td>0.409</td>\n","      <td>0.92975</td>\n","      <td>1.295</td>\n","      <td>1.8425</td>\n","      <td>4.84</td>\n","    </tr>\n","    <tr>\n","      <th>LC</th>\n","      <td>2016.0</td>\n","      <td>1.351375</td>\n","      <td>0.645456</td>\n","      <td>0.402</td>\n","      <td>0.87775</td>\n","      <td>1.200</td>\n","      <td>1.6900</td>\n","      <td>4.59</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["       count      mean       std    min      25%    50%     75%   max\n","conf                                                                 \n","HC    1972.0  1.454450  0.706737  0.409  0.92975  1.295  1.8425  4.84\n","LC    2016.0  1.351375  0.645456  0.402  0.87775  1.200  1.6900  4.59"]},"execution_count":15,"metadata":{},"output_type":"execute_result"}],"source":["# 调用pandas对象自带的方法，将数据按照‘conf’条件进行分组，选择每组的反应时间，展示描述性统计的内容\n","data.groupby('conf').rt.describe() "]},{"cell_type":"code","execution_count":16,"metadata":{"collapsed":false,"id":"ECD0D587538146E590F20C9CD4C83390","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["<AxesSubplot:ylabel='Density'>"]},"execution_count":16,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/ECD0D587538146E590F20C9CD4C83390/rkyxu3gf20.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# 调用pandas对象自带的画图功能，选择每组的反应时间，绘制概率密度函数图\n","data.rt.plot.density() "]},{"cell_type":"code","execution_count":17,"metadata":{"collapsed":false,"id":"A535064986504AAC999A0168FCE6DD6C","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["conf\n","HC    AxesSubplot(0.125,0.125;0.775x0.755)\n","LC    AxesSubplot(0.125,0.125;0.775x0.755)\n","Name: rt, dtype: object"]},"execution_count":17,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/A535064986504AAC999A0168FCE6DD6C/rkyxu3ptlo.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# 调用pandas对象自带的画图功能，按照‘conf’进行分组，选择每组的反应时间，绘制概率密度函数图\n","data.groupby(['conf']).rt.plot.density()"]},{"cell_type":"code","execution_count":18,"metadata":{"collapsed":false,"id":"12A3713F9D63424D8E3F444703C04B24","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["# 将‘conf’条件进行0， 1编码，HC转化为1，LC转化为0，以便在pymc3中进行计算\n","data.conf = data.conf.map({'HC':1,'LC':0})"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"F927A94BBCD540508F38A653E8CA8EA9","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### Model1: Normal\n","\n","首先，我们假定反应时是呈正态分布的，正态分布的均值是由高一致性条件和低一致性条件决定的。\n","\n","\n","![Image Name](https://docs.pymc.io/en/v3/_images/continuous-18.png)\n","\n","图片来源：https://docs.pymc.io/en/v3/api/distributions/continuous.html#pymc3.distributions.continuous.Normal"]},{"cell_type":"code","execution_count":19,"metadata":{"collapsed":false,"id":"96E1B6A51A5E4995AF3E0D391A1E6978","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with pm.Model() as NormalModel:\n","    # 先验分布: alpha, beta, sigma这三个参数是随机变量\n","    sigma = pm.HalfNormal('sigma', sd=1)\n","    alpha = pm.Normal('alpha', mu=0, sd=1)\n","    beta = pm.Normal('beta', mu=0, sd=1)\n","    # 自变量conf是之前已经载入的数据\n","    x = pm.Data(\"x\", data['conf'])\n","    # 正态分布均值是确定性随机变量，这个变量的值完全由右端值确定\n","    mu = pm.Deterministic(\"mu\", alpha + beta*x) \n","    # Y的观测值，这是一个特殊的观测随机变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n","    y_obs = pm.Normal('y_obs',mu=mu,sd=sigma,observed=data['rt'] )"]},{"cell_type":"code","execution_count":20,"metadata":{"collapsed":false,"id":"0BDF76D3D68F4C78A21B3361B98C364E","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/0BDF76D3D68F4C78A21B3361B98C364E/rkyxu5eeme.svg\">"],"text/plain":["<graphviz.graphs.Digraph at 0x7f41c607f5d0>"]},"execution_count":20,"metadata":{},"output_type":"execute_result"}],"source":["# 展示模型结构\n","pm.model_to_graphviz(NormalModel)"]},{"cell_type":"code","execution_count":21,"metadata":{"collapsed":false,"id":"EC6B63389E8D4AFB8CA5450DE5D8D4A9","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Multiprocess sampling (2 chains in 2 jobs)\n","NUTS: [beta, alpha, sigma]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"name":"stderr","output_type":"stream","text":["Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 7 seconds.\n"]}],"source":["with NormalModel:\n","    # 使用mcmc方法进行采样，draws为采样次数，tune为调整采样策略的次数，这些次数将在采样结束后被丢弃，\n","    # target_accept为接受率， return_inferencedata=True为该函数返回的对象是arviz.InnferenceData对象\n","    # chains为我们采样的链数，cores为我们的调用的cpu数，多个链可以在多个cpu中并行计算，我们在和鲸中调用的cpu数为2\n","    trace = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)"]},{"cell_type":"code","execution_count":22,"metadata":{"collapsed":false,"id":"72F6BC0AFD5E4BFBBE90BE84126D8A3B","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["array([[<AxesSubplot:title={'center':'alpha'}>,\n","        <AxesSubplot:title={'center':'alpha'}>],\n","       [<AxesSubplot:title={'center':'beta'}>,\n","        <AxesSubplot:title={'center':'beta'}>],\n","       [<AxesSubplot:title={'center':'sigma'}>,\n","        <AxesSubplot:title={'center':'sigma'}>],\n","       [<AxesSubplot:title={'center':'mu'}>,\n","        <AxesSubplot:title={'center':'mu'}>]], dtype=object)"]},"execution_count":22,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/72F6BC0AFD5E4BFBBE90BE84126D8A3B/rkyy3ipr99.png\">"],"text/plain":["<Figure size 864x576 with 8 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# 绘制各参数的采样情况\n","az.plot_trace(trace,var_names=['alpha','beta','sigma','mu'])"]},{"cell_type":"code","execution_count":23,"metadata":{"collapsed":false,"id":"0950B0FA00DB48DD98BE536B08BA8DBB","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>mean</th>\n","      <th>sd</th>\n","      <th>hdi_3%</th>\n","      <th>hdi_97%</th>\n","      <th>mcse_mean</th>\n","      <th>mcse_sd</th>\n","      <th>ess_bulk</th>\n","      <th>ess_tail</th>\n","      <th>r_hat</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>alpha</th>\n","      <td>1.351</td>\n","      <td>0.015</td>\n","      <td>1.323</td>\n","      <td>1.380</td>\n","      <td>0.000</td>\n","      <td>0.0</td>\n","      <td>1749.0</td>\n","      <td>2413.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>beta</th>\n","      <td>0.104</td>\n","      <td>0.021</td>\n","      <td>0.062</td>\n","      <td>0.142</td>\n","      <td>0.001</td>\n","      <td>0.0</td>\n","      <td>1469.0</td>\n","      <td>1686.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>sigma</th>\n","      <td>0.677</td>\n","      <td>0.008</td>\n","      <td>0.662</td>\n","      <td>0.691</td>\n","      <td>0.000</td>\n","      <td>0.0</td>\n","      <td>2621.0</td>\n","      <td>2656.0</td>\n","      <td>1.0</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["        mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n","alpha  1.351  0.015   1.323    1.380      0.000      0.0    1749.0    2413.0   \n","beta   0.104  0.021   0.062    0.142      0.001      0.0    1469.0    1686.0   \n","sigma  0.677  0.008   0.662    0.691      0.000      0.0    2621.0    2656.0   \n","\n","       r_hat  \n","alpha    1.0  \n","beta     1.0  \n","sigma    1.0  "]},"execution_count":23,"metadata":{},"output_type":"execute_result"}],"source":["# 参数的统计值\n","az.summary(trace,var_names=['alpha','beta','sigma'])"]},{"cell_type":"code","execution_count":24,"metadata":{"collapsed":false,"id":"59EAD4BE19CC48CD86F069F10D8823BD","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"data":{"text/html":["\n","    <div>\n","      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n","      100.00% [4000/4000 00:06&lt;00:00]\n","    </div>\n","    "],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"update_display_data"},{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n","  FutureWarning,\n"]}],"source":["with NormalModel:\n","     # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n","    ppc_y = pm.sample_posterior_predictive(trace.posterior) \n","#将ppc_y转化为InferenceData对象合并到trace中\n","az.concat(trace, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)"]},{"cell_type":"code","execution_count":25,"metadata":{"collapsed":false,"id":"735E869AB9864E0887F4304F6CE4C614","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["<AxesSubplot:xlabel='y_obs'>"]},"execution_count":25,"metadata":{},"output_type":"execute_result"},{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/IPython/core/events.py:89: UserWarning: Creating legend with loc=\"best\" can be slow with large amounts of data.\n","  func(*args, **kwargs)\n"]},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/735E869AB9864E0887F4304F6CE4C614/rkyy4ne1fn.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# 绘制后验预测分布\n","az.plot_ppc(trace)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"D7A33179EAC94F968EBEDCE18B9F7B54","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### Model2: Log-normal\n","\n","然后，我们假定反应时是呈log正态分布的，LogNormal分布的均值是由高一致性条件和低一致性条件决定的。\n","\n","\n","![Image Name](https://docs.pymc.io/en/v3/_images/continuous-16.png)\n","\n","\n","\n","图片来源：https://docs.pymc.io/en/v3/api/distributions/continuous.html#pymc3.distributions.continuous.LogitNormal"]},{"cell_type":"code","execution_count":26,"metadata":{"collapsed":false,"id":"55D72B3FF6274317829D74EC8C4C5F99","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with pm.Model() as LogNormal:\n","    # 先验分布: alpha, beta, sigma这三个参数是随机变量\n","    sigma = pm.HalfNormal('sigma', sd=1)\n","    alpha = pm.Normal('alpha', mu=0, sd=1)\n","    beta = pm.Normal('beta', mu=0, sd=1)\n","    # 自变量conf是之前已经载入的数据\n","    x = pm.Data(\"x\", data['conf'])\n","    # 正态分布均值是确定性随机变量，这个变量的值完全由右端值确定\n","    mu = pm.Deterministic(\"mu\", alpha + beta*x) \n","    # Y的观测值，这是一个特殊的观测随机变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n","    y_obs = pm.Lognormal('y_obs',mu=mu,sd=sigma,observed=data['rt'] )"]},{"cell_type":"code","execution_count":27,"metadata":{"collapsed":false,"id":"30736F5A3C6B499282B9CC9D5426D679","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/30736F5A3C6B499282B9CC9D5426D679/rkyy4rgv8l.svg\">"],"text/plain":["<graphviz.graphs.Digraph at 0x7f41a56a4f50>"]},"execution_count":27,"metadata":{},"output_type":"execute_result"}],"source":["# 展示模型结构\n","pm.model_to_graphviz(LogNormal)"]},{"cell_type":"code","execution_count":28,"metadata":{"collapsed":false,"id":"4461A25C102D41CB87894A78FACFE2F2","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Multiprocess sampling (2 chains in 2 jobs)\n","NUTS: [beta, alpha, sigma]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"name":"stderr","output_type":"stream","text":["Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 6 seconds.\n"]}],"source":["with LogNormal:\n","    # 使用mcmc方法进行采样，draws为采样次数，tune为调整采样策略的次数，这些次数将在采样结束后被丢弃，\n","    # target_accept为接受率， return_inferencedata=True为该函数返回的对象是arviz.InnferenceData对象\n","    # chains为我们采样的链数，cores为我们的调用的cpu数，多个链可以在多个cpu中并行计算，我们在和鲸中调用的cpu数为2\n","    trace2 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)"]},{"cell_type":"code","execution_count":29,"metadata":{"collapsed":false,"id":"C37FA075D48842BB8EF5C3ED450B8954","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["array([[<AxesSubplot:title={'center':'alpha'}>,\n","        <AxesSubplot:title={'center':'alpha'}>],\n","       [<AxesSubplot:title={'center':'beta'}>,\n","        <AxesSubplot:title={'center':'beta'}>],\n","       [<AxesSubplot:title={'center':'sigma'}>,\n","        <AxesSubplot:title={'center':'sigma'}>],\n","       [<AxesSubplot:title={'center':'mu'}>,\n","        <AxesSubplot:title={'center':'mu'}>]], dtype=object)"]},"execution_count":29,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/C37FA075D48842BB8EF5C3ED450B8954/rkyyeb1nx8.png\">"],"text/plain":["<Figure size 864x576 with 8 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# 绘制各参数的采样情况\n","az.plot_trace(trace2,var_names=['alpha','beta','sigma','mu'])"]},{"cell_type":"code","execution_count":30,"metadata":{"collapsed":false,"id":"6DFF8D88AFE64DEE8382C5BA63734C7E","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>mean</th>\n","      <th>sd</th>\n","      <th>hdi_3%</th>\n","      <th>hdi_97%</th>\n","      <th>mcse_mean</th>\n","      <th>mcse_sd</th>\n","      <th>ess_bulk</th>\n","      <th>ess_tail</th>\n","      <th>r_hat</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>alpha</th>\n","      <td>0.197</td>\n","      <td>0.010</td>\n","      <td>0.178</td>\n","      <td>0.217</td>\n","      <td>0.0</td>\n","      <td>0.0</td>\n","      <td>2051.0</td>\n","      <td>2307.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>beta</th>\n","      <td>0.066</td>\n","      <td>0.015</td>\n","      <td>0.040</td>\n","      <td>0.095</td>\n","      <td>0.0</td>\n","      <td>0.0</td>\n","      <td>2124.0</td>\n","      <td>2102.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>sigma</th>\n","      <td>0.465</td>\n","      <td>0.005</td>\n","      <td>0.455</td>\n","      <td>0.475</td>\n","      <td>0.0</td>\n","      <td>0.0</td>\n","      <td>2211.0</td>\n","      <td>2175.0</td>\n","      <td>1.0</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["        mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n","alpha  0.197  0.010   0.178    0.217        0.0      0.0    2051.0    2307.0   \n","beta   0.066  0.015   0.040    0.095        0.0      0.0    2124.0    2102.0   \n","sigma  0.465  0.005   0.455    0.475        0.0      0.0    2211.0    2175.0   \n","\n","       r_hat  \n","alpha    1.0  \n","beta     1.0  \n","sigma    1.0  "]},"execution_count":30,"metadata":{},"output_type":"execute_result"}],"source":["# 参数的统计值\n","az.summary(trace2,var_names=['alpha','beta','sigma'])"]},{"cell_type":"code","execution_count":31,"metadata":{"collapsed":false,"id":"2E69A3EAF68846CD95EC6C15E5C4BFB8","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"data":{"text/html":["\n","    <div>\n","      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n","      100.00% [4000/4000 00:05&lt;00:00]\n","    </div>\n","    "],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"update_display_data"},{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n","  FutureWarning,\n"]}],"source":["with LogNormal:\n","    # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n","    ppc_y = pm.sample_posterior_predictive(trace2.posterior) \n","#将ppc_y转化为InferenceData对象合并到trace中\n","az.concat(trace2, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)"]},{"cell_type":"code","execution_count":32,"metadata":{"collapsed":false,"id":"84AA12B197AE455D8BF40FD167BF84AA","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["<AxesSubplot:xlabel='y_obs'>"]},"execution_count":32,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/84AA12B197AE455D8BF40FD167BF84AA/rkyyfcyzbu.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# 绘制后验预测分布\n","az.plot_ppc(trace2)"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"676785E87AB84CCA8713AB6F1AAC6255","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["#### Model3: Gamma\n","\n","最后，我们假定反应时是呈gamma分布的，gamma分布的alpha是由高一致性条件和低一致性条件决定的。\n","\n","\n","![Image Name](https://docs.pymc.io/en/v3/_images/continuous-6.png)\n","\n","图片来源：https://docs.pymc.io/en/v3/api/distributions/continuous.html#pymc3.distributions.continuous.Gamma\n"]},{"cell_type":"code","execution_count":33,"metadata":{"collapsed":false,"id":"0120ACDED2BA4D2F8631B14B2A512624","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":["with pm.Model() as Gamma:\n","    # 先验分布:intercept, theta, beta这三个参数是随机变量\n","    intercept = pm.HalfNormal('intercept', sd=1)\n","    theta = pm.HalfNormal('theta', sd=1)\n","    beta = pm.HalfNormal('beta', sd=1)   \n","    # 自变量conf是之前已经载入的数据\n","    x = pm.Data(\"x\", data['conf'])\n","    # 参数alpha是确定性随机变量，这个变量的值完全由右端值确定\n","    alpha = pm.Deterministic(\"alpha\",  intercept+ theta*x) \n","    # Y的观测值，这是一个特殊的观测随机变量，表示模型数据的可能性。也可以表示模型的似然，通过 observed 参数来告诉这个变量其值是已经被观测到了的，不会被拟合算法改变\n","    y_obs = pm.Gamma('y_obs',alpha=alpha,beta=beta,observed=data['rt'] )"]},{"cell_type":"code","execution_count":34,"metadata":{"collapsed":false,"id":"BE80C3787B36411C8A0B2893446C57B7","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/BE80C3787B36411C8A0B2893446C57B7/rkyyfor753.svg\">"],"text/plain":["<graphviz.graphs.Digraph at 0x7f41695b42d0>"]},"execution_count":34,"metadata":{},"output_type":"execute_result"}],"source":["# 展示模型结构\n","pm.model_to_graphviz(Gamma)"]},{"cell_type":"code","execution_count":35,"metadata":{"collapsed":false,"id":"595F516FDCAB4A3C821F09CE6B8FE28D","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Multiprocess sampling (2 chains in 2 jobs)\n","NUTS: [beta, theta, intercept]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"name":"stderr","output_type":"stream","text":["Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 24 seconds.\n"]}],"source":["with Gamma:\n","    # 使用mcmc方法进行采样，draws为采样次数，tune为调整采样策略的次数，这些次数将在采样结束后被丢弃，\n","    # target_accept为接受率， return_inferencedata=True为该函数返回的对象是arviz.InnferenceData对象\n","    # chains为我们采样的链数，cores为我们的调用的cpu数，多个链可以在多个cpu中并行计算，我们在和鲸中调用的cpu数为2\n","    trace3 = pm.sample(draws = 2000, tune=1000, target_accept=0.9,chains=2, cores= 2,return_inferencedata=True)"]},{"cell_type":"code","execution_count":36,"metadata":{"collapsed":false,"id":"7D89712978D347B98A92A0F9ACF7D28D","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["array([[<AxesSubplot:title={'center':'beta'}>,\n","        <AxesSubplot:title={'center':'beta'}>],\n","       [<AxesSubplot:title={'center':'intercept'}>,\n","        <AxesSubplot:title={'center':'intercept'}>],\n","       [<AxesSubplot:title={'center':'theta'}>,\n","        <AxesSubplot:title={'center':'theta'}>]], dtype=object)"]},"execution_count":36,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/7D89712978D347B98A92A0F9ACF7D28D/rkyyh4r2l.png\">"],"text/plain":["<Figure size 864x432 with 6 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# 绘制各参数的采样情况\n","az.plot_trace(trace3,var_names=['beta','intercept','theta'])"]},{"cell_type":"code","execution_count":37,"metadata":{"collapsed":false,"id":"0F6385F002B14CA9805B85AB20D4F94F","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>mean</th>\n","      <th>sd</th>\n","      <th>hdi_3%</th>\n","      <th>hdi_97%</th>\n","      <th>mcse_mean</th>\n","      <th>mcse_sd</th>\n","      <th>ess_bulk</th>\n","      <th>ess_tail</th>\n","      <th>r_hat</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>beta</th>\n","      <td>3.358</td>\n","      <td>0.076</td>\n","      <td>3.220</td>\n","      <td>3.501</td>\n","      <td>0.002</td>\n","      <td>0.002</td>\n","      <td>1230.0</td>\n","      <td>1515.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>intercept</th>\n","      <td>4.569</td>\n","      <td>0.103</td>\n","      <td>4.376</td>\n","      <td>4.755</td>\n","      <td>0.003</td>\n","      <td>0.002</td>\n","      <td>1192.0</td>\n","      <td>1627.0</td>\n","      <td>1.0</td>\n","    </tr>\n","    <tr>\n","      <th>theta</th>\n","      <td>0.286</td>\n","      <td>0.066</td>\n","      <td>0.165</td>\n","      <td>0.408</td>\n","      <td>0.002</td>\n","      <td>0.001</td>\n","      <td>1864.0</td>\n","      <td>1289.0</td>\n","      <td>1.0</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["            mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\n","beta       3.358  0.076   3.220    3.501      0.002    0.002    1230.0   \n","intercept  4.569  0.103   4.376    4.755      0.003    0.002    1192.0   \n","theta      0.286  0.066   0.165    0.408      0.002    0.001    1864.0   \n","\n","           ess_tail  r_hat  \n","beta         1515.0    1.0  \n","intercept    1627.0    1.0  \n","theta        1289.0    1.0  "]},"execution_count":37,"metadata":{},"output_type":"execute_result"}],"source":["# 参数的统计值\n","az.summary(trace3,var_names=['beta','intercept','theta'])"]},{"cell_type":"code","execution_count":38,"metadata":{"collapsed":false,"id":"AE83F205FC3C45C782E4767118C39CA4","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"data":{"text/html":["\n","    <div>\n","      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n","      100.00% [4000/4000 00:07&lt;00:00]\n","    </div>\n","    "],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"update_display_data"},{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/arviz/data/io_pymc3.py:100: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.\n","  FutureWarning,\n"]}],"source":["with Gamma:\n","    # pm.sample_posterior_predictive()利用trace.posterior的后验分布计算后验预测分布\n","    ppc_y = pm.sample_posterior_predictive(trace3.posterior) \n","#将ppc_y转化为InferenceData对象合并到trace中\n","az.concat(trace3, az.from_pymc3(posterior_predictive=ppc_y), inplace=True)"]},{"cell_type":"code","execution_count":39,"metadata":{"collapsed":false,"id":"4C65964DA6FC4648819B8327252C3BC8","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"data":{"text/plain":["<AxesSubplot:xlabel='y_obs'>"]},"execution_count":39,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/4C65964DA6FC4648819B8327252C3BC8/rkyyi9cbig.png\">"],"text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["# 绘制后验预测分布\n","az.plot_ppc(trace3)"]},{"cell_type":"code","execution_count":40,"metadata":{"collapsed":false,"id":"42C27ED05D8F4F698BBC58DE8F6AE318","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[{"name":"stderr","output_type":"stream","text":["/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n","  \"The default method used to estimate the weights for each model,\"\n"]},{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>rank</th>\n","      <th>loo</th>\n","      <th>p_loo</th>\n","      <th>d_loo</th>\n","      <th>weight</th>\n","      <th>se</th>\n","      <th>dse</th>\n","      <th>warning</th>\n","      <th>loo_scale</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>log-nomal</th>\n","      <td>0</td>\n","      <td>-3520.882658</td>\n","      <td>2.803835</td>\n","      <td>0.000000</td>\n","      <td>1.000000e+00</td>\n","      <td>49.737714</td>\n","      <td>0.000000</td>\n","      <td>False</td>\n","      <td>log</td>\n","    </tr>\n","    <tr>\n","      <th>gamma</th>\n","      <td>1</td>\n","      <td>-3596.278905</td>\n","      <td>2.774815</td>\n","      <td>75.396247</td>\n","      <td>0.000000e+00</td>\n","      <td>48.999506</td>\n","      <td>8.475056</td>\n","      <td>False</td>\n","      <td>log</td>\n","    </tr>\n","    <tr>\n","      <th>normal</th>\n","      <td>2</td>\n","      <td>-4102.280133</td>\n","      <td>3.805405</td>\n","      <td>581.397475</td>\n","      <td>1.440962e-09</td>\n","      <td>58.326281</td>\n","      <td>29.497056</td>\n","      <td>False</td>\n","      <td>log</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["           rank          loo     p_loo       d_loo        weight         se  \\\n","log-nomal     0 -3520.882658  2.803835    0.000000  1.000000e+00  49.737714   \n","gamma         1 -3596.278905  2.774815   75.396247  0.000000e+00  48.999506   \n","normal        2 -4102.280133  3.805405  581.397475  1.440962e-09  58.326281   \n","\n","                 dse  warning loo_scale  \n","log-nomal   0.000000    False       log  \n","gamma       8.475056    False       log  \n","normal     29.497056    False       log  "]},"execution_count":40,"metadata":{},"output_type":"execute_result"}],"source":["# 将三个模型的采样结果进行比较\n","compare_dict = {\"normal\": trace, \"log-nomal\": trace2, \"gamma\": trace3}\n","# 选择loo方法进行比较\n","comp = az.compare(compare_dict, ic='loo')\n","comp"]},{"cell_type":"markdown","metadata":{"collapsed":false,"id":"7F6A2002F94D44FC9B299B76C86EE6B9","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":["- 第一列为索引，它列出了传递给 az.compare(.) 的模型名称。\n","\n","- rank 列：按照预测精度做的排名，值从0依次到模型总数，其中0代表最高精度。\n","\n","- loo 列：各模型 ELPD 值的列表，总是按照 ELPD 值从最好到最差排序。\n","\n","- p_loo 列：惩罚项的值列表，可以将其粗略地视为有效参数数量的估计值（但不要太认真）。此值可能低于具有更多结构的模型（如分层模型）中的实际参数数量，或者高于那些预测能力非常弱或严重错误指定的模型的实际参数数量。\n","\n","- d_loo 列：每个模型与排名第一的模型之间的 LOO 相对差。因此第一个模型始终取值为0  。\n","\n","- weight 列：分配给每个模型的权重。权重可以粗略地解释为在指定数据的条件下，是（参与比较的各模型中）该模型的概率。\n","\n","- se 列：ELPD 的标准误差。\n","\n","- dse 列：ELPD 相对差的标准误差。 dse 与 se 不一定相同，因为 ELPD 的不确定性在模型之间可能存在相关性。排名第一的模型 dse 值始终为0 。\n","\n","- warning 列：如果为True，表示这是一个警告，LOO 的近似估计不可靠。\n","\n","- loo_scale 列：估计值所用的尺度。默认为对数尺度。其他选项还包括：离差值尺度，即对数分值乘以-2，这会颠倒排序，较低的 ELPD 会更好；负对数尺度，即对数分值乘以-1，与离差值尺度一样，值越低越好。"]},{"cell_type":"code","execution_count":null,"metadata":{"collapsed":false,"id":"0D3923A1486941098AC9741BFEE5BFDB","jupyter":{},"notebookId":"6368be22d28b18529a659034","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"outputs":[],"source":[]}],"metadata":{"kernelspec":{"display_name":"Python 3.9.6 64-bit","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.9.6"},"vscode":{"interpreter":{"hash":"31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"}}},"nbformat":4,"nbformat_minor":0}
